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We theoretically study current dynamics of graphene nanoribbons subject to bias dc and ac driven fields. 
We showed that graphene nanoribbons exhibit negative high-harmonic differential conductivity. Negative 
differential conductivity appears when bias electric filed is in the neighborhood of applied ac filed amplitude. 
We also observe both even and odd high-harmonic negative differential conductivity at wave mixing of two 
commensurate frequencies. The even harmonics are more pronounced than the odd harmonics. A possible 
use of the present method for generating terahertz frequencies at even harmonics in graphene is suggested. 
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I. INTRODUCTION 

Graphene has continued to surprise scientists since its 
discovery in 2004 by Geim and his tearrP} Theoretically, 
the carrier transport properties are fantastic. Especially, 
its high carrier mobility of 44000cm 2 V r_1 s _1 '2l. But at- 
tempts to utilize these astonishing properties in graphene 
devices are posing some difficulties. The limitation is 
probably due to several factors including; lack of band 
gap in graphene sheets, edge defects, disorder, among 
others. To overcome some of these obstacles, particu- 
larly the lack of band gap, the dimension of graphene 
sheet has to be reduced. After all, new physics (quanti- 
zation) emerge when dimensions of materials are reduced. 
An infinite 2D graphene could become a one dimensional 
plus quantization along one other direction, the con- 
sequence is the opening of a gap in the its dispersion 
spectrum. The resulting material from the 2D infinite 
sheet form a graphene nanoribbon (GNR). Depending 
on the nature of the ribbon edges, one gets two symme- 
try groups; Armchair Graphene Nanoribbon (aGNR) and 
Zigzag Graphene Nanoribbon (zGNR). Electron dynam- 
ics of both aGNR and zGNR have different properties, 
mostly due to the berry phase and pseudo spirP. Edge 
states have significant contribution to graphene proper- 
ties, because in a nanometer size ribbon, massless Dirac 
fermions can reach the edges within a femto-second be- 
fore encountering any other lattice effects, like electron- 
electron interaction, electron-phonon interaction, etc. 

In this paper, we study the phenomenon of negative 
differential conductivity (NDC) in GNRs. In conven- 
tional semiconductor devices, a negative differential con- 
ductive behavior is known to offer great potential for high 
frequency applications as Bloch oscillators, frequency 
multipliers, and fast switching devices. For this reason, 
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the NDC effect has been greatly explored and discussed in 
several graphene nanostructures, particularly irP. NDC 
can also be observed in other graphene allotropes; car- 
bon nanotubes (CNT^l The unique energy spectrum of 
holes and electrons in GNRs, especially its narrow gap- 
less nature leads to nontrivial features such as Negative 
Differential Conductivity (NDC) in the THz regim<P. 

In fact, we must emphasize that the phenomenon of 
NDC and Bloch oscillations in a material is a possible 
signature for THz generations in the materia l 7 * 11 !, since 
NDC occurs in the THz spectral range. Most of the 
methods of producing THz frequencies are experimen- 
tal and only very few analytical (without computer nu- 
merics) approaches are known. Though Green function 
techniques have also been employed in some cases. Mo- 
tivated by the fact that a rigorous analytical approach is 
necessary for studying NDC in GNRs, we adopt a semi- 
classical method used in reference d 5 * 12 ! for armchair and 
zigzag CNTs. We predict that GNRs should also repro- 
duce similar results as in reference^ because both CNTs 
and GNRs have almost the same full tight binding (TB) 
nonlinear complex band structure. 

The rest of this paper is organized as follows; In Ssec- 
tion [TTJ we derived the current density of aGNR and 
zGNR and imposed certain conditions to reduce the 
equations to forms appropriate for our models. The re- 
sults obtained in Section |III| are plotted and discussed 
in Section |IV| The paper finally concludes in Section [V] 
where some recommendations for future applications are 
made. 



II. THE THEORY 

As it is usually done in semi-classical treatment of 
quantum systems, we assume that the dynamics of the 
free 7r-electrons in graphene satisfies the time-dependent 
Boltzmann transport equation (BTE) in zero magnetic 
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field. That is, 

df(k,t) , eE(t)df(k,t) 



at 



dk 



T[f(k,t) - f (k)]. (1) 



We are also assuming relaxation time approximation 
(RTA) and a spatial uniform graphene nanoribbon. Also, 
the inverse of the relaxation time T is momentum inde- 
pendent. For the case of energy varying T, se e 3 . In EqjlJ 
fo(k) and f(t, k) are the equilibrium and non-equilibrium 
Fermi electron distribution functions, respectively, e is 
the electronic charge, k is the electron wave vector and h 
is the reduced plank constant. We consider an external 
applied field 



(2) 



3=0 



as a superposition of n harmonic waves polarized along 
one direction with the angular frequency uj. The phase 
difference between the (j + l)th and jth wave being 
oy+x — Q.j = a is arbitrary, j is an integer. Ej are the 
field amplitudes. We require that ujq = ccq = 0. In 
the following section we will look at current density for 
GNRs. 



A. Armchair and zigzag nanoribbon band structures 

The energy band structure of aGNR and zGNR is char- 
acterized by t hree parameters; band index A, phase 9 and 
wave vector aP^. For aGNR 



S x (k,9) = A 7ov /l + 4cos 2 {sA9) + 4cos(sA6)cos(kl), 



and for zGNR 



(3) 



S x (k,9) = A 7ov /l + Acos 2 {kV) + 4cos{sA6)cos(kl'). 

(4) 

A = ±1. (+) for conduction band and (-) for valence 
band. I = v3o/2 and I' — a/2, a is the lattice spacing 
with value 0.246nm, 70 ~ S.OeV is the overlap integral 
and 9 is the phase perpendicular to the quasi-momentum 
hk. The 1BZ of aGNR is bounded by kl = [-7r/2,7r/2] 
and the zGNR is kl' — [0,tt]. k is parallel to the edge 
and has translational symmetry along this direction. For 
aGNR, the transverse wave vector (phase) is quantized 
according to the rule^, 9 S = sA9 with A9 = an d 
s = 1,2,... , j\f. Unlike aGNR, the nature of transverse 
wave vector quantization is complicated in zGNR, de- 
pending on both k and 9 as A9 S = (7rj+A(fc, 9))/(N+l). 
However, for simplicity we assume A is constant, say 7r/2, 

so that A6 S = ^ 2 js+1* ■ Except this little subtlety for 
zGNRs, all that will be discussed in the following for 
aGNR are equally applied to the zGNR. 

Now, employed translational invariance of the 
graphene ribbon in the reciprocal space and expand in 



Fourier series functions /, /o and £ along the edge hav- 
ing the periodicity in k. i.e, 



/o(M) = 



irkl 



f(k,9,t)=Y / fr(0y rkl $r(t), 



£(k,6) 



7o2j£ r (0)e 



irkl 



(5) 



(6) 



(7) 



The Fourier coefficient f r is expressed as f r (9) = 
UsA95{9 s - sA9) with 
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Irs 



and 



■ksAO 



I 



tt/21 



dkf (k,9)e 



— irkl 



2tt7o 



tt/21 



t/21 



dk£{k)e 



-irkl 



frs = I' 



£ r = £*_ 



(8) 



(9) 



s in Eq.Q counts the number of dimers j\f in GNRs. 
The factor $ r in Eq.Q is a central point in this paper 
and so has to be determined, r is an integer and not 
equal to zero. We consider a classical limit in which en- 
ergy levels could be excited due to thermal fluctuations, 
i.e AS « KbT « Sc- This condition is also neces- 
sary for large enough field, so that charge carriers can 
escape low energy scattering^. The energy level spacing 
AS = ■nW^ol/A, Sc is the charging energy, Kb is the 
Boltzmann constant, T is the lattice temperature and 
W is the graphene width. In what follows next, we will 
find the form of To do this Eqs.([5]), (|6|, are 

substituted in Eq.Q to yield 



dt 



[T + irfi(t)]$ r (f) 



0, 



(10) 



where Q(t) = fl + f Eje lu >J t+a i is the modulation 

degree of anharmonicity in electron motion. The solution 
of eq.(10) is a straight forward one, using the boundary 
conditions t = 0, $ r (0) = 1, one has 

T f die rt+J £?=i ft ' ei " J ' t+Q3+i/3ot 
®r(t) = — i — tz , (11) 

f3 j = erlEj /Hujj and /3q = Qo — erlE$. One can intro- 
duce product notation in eq.(ll) as 



ne<^ r (t) = r f[ n \e- rt - i ^ cos( - iu 'i t+a '^- iflot 

x / di e rt + i f } j cos ( i ^j t + a j)+ i l 3 o t (\2) 
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Zm$ r (t) = T II II [e -rt- ^^* iri K*+ a! i)-^o' 
i'=ij=i 

X / ( j i f e rt + i 0i sin ( iu 'i t + a i)+ i l 3 a t ^3) 



Eqs.(12|, (13) are connected with the well known Bessel 



functions via Jacobi- Anger expansion 



giving 



with 



3(t) =i^2jo,r®r(t) + C.C 



(19) 



30,r 



Jo. 



~JO,-r- 



±im6 



m— — oo 



cosQ 



E i m Jm(fr)t 



Jm{f3) is the m th order Bessel function. Using the ex- 
pansions above, Eq.(13) becomes 



oo oo n n 



Tij/— — oo rrij — — ooj' — lj — l 
x g— rt— iPjisintiu^t+ciji)— iftjt 

x / dtJ m3 {P ] )e Tt+l ^ sm{lu ' t+ai)+tf)ot (14) 



Applying the integration in the preceding equation and 
letting j = f = 1, 2, 3, . . . and rrij, rij — ±1, ±2, ±3, . . . 



oo n 



rij , i/j — — oo j — 1 



where Vj = — rrij . 



OJj t-\-ilSj Otj 

+ zr(/3 + " 
(15) 



Substitute $ r in Eq.(|19| to get 



oo oo n 

j(t) = »X>,r E IT J ", (M J n } -^ (&) 
r—1 rij, Vj— — oo j — 1 

x 1 ZTTfl 1 V + c - c - ( 2 °) 

Note the r dependence of /3j, (3q and the summation 
over the index. Using the formalism by Litvinov and 
ManassorfS Eq.(20l can be put in a taylor-like expan- 
sion of Ej. i.e 

m = j dc + \ e E i E ^w^-* +«+-, (2i) 



where 



oo n 



Jdc 



E '- E iI<^) T 

r—1 rij—~coj — l 



i + (3()T + Tl,jUJjT 

+ [r(/3 + n J ^)] 2 



is the differential dc conductivity (for tA, = 0), and 



(22) 



oo n 



<r„, 2E./,,, E II JnMJ ^~ M) 

r—1 rij— — ocj — l 3 

l + [r(/3 + n^)] 2 

is the large-signal dynamic nonlinear conductivity at Vj 
harmonic with drive frequency ujj . 



B. Sheet current density 

The sheet current density of graphene can be deter- 
mined from the relation 



III. NEGATIVE DIFFERENTIAL CONDUCTIVITY 



A. Pure dc limit 



■ 9s9v v-v , 

3it) = — 2_^ev k f k 

k 



(16) 



The sheet area A = WL, with W the width. g s , g v are 
the spin and valley degeneracies respectively. For the 
aGNR, 



m 



g s g v e 

4tt 2 



E J dkv(k,e s )f(k,6 s ,$ r ). (17) 



The velocity of Dirac fermions in graphene is defined as 
v(k) = d£/hdk. In terms of the Fourier coefficients, 



rkl 



(18) 



To see immediately that Eq.(22| demonstrates NDC, 



we consider a pure dc limit where Uj — > 0. The 
Bessel functions except the n = term will vanish. 
The real part of the differential conductivity cr(0) = 
lim Uj ._ i .o dj/dE becomes 



(24) 



so that if /3q > t -1 , the differential conductivity is nega- 
tive and NDC is manifest in GNRs. 

Electron dynamics may be come more complicated in 
the presence of high-frequency components in addition 
to the static electric fields. High negative differential 
conductivity thus may result in GNRs if an external drive 
force is applied. 
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B. Monoharmonics 



If one component of an ac filed in Eq.|l]) is applied, 
then n = 1 and Eq.(22| simplifies to 



f3 r + nuiT 
+ (/3 t + ulot) 2 ' 



(25) 



after dropping the subscripts on n. Here, it not clear 
immediately how NDC can be seen. To observe it, we 
plot j versus Eq for lot « 1 is shown in Fig[l](left) for 
armchair ribbon and in Fig[l] (right) for zigzag ribbon, 



with fix. 2 = erlE\^jhbj\^. Where I = \fZaj1 for arm- 
chair and V = a/2 for TA'gL&'g graphene nanoribbons re- 
spectively. Finally, the current density becomes 



eh 2 



E 2 cosa 



for 



TLLOT 



+ (fi T + TILOt) 2 



xJ n (rfix) {Jn-^rpx) - J„ + „(rfix)} , (29) 

which reduces to the monoharmonic case when n = 1. 

The nature of the NDC is observed for a simultaneously 
varying harmonic field and phase difference in a three 
dimensional plot shown in Fig(3] 



C. Biharmonics 

One can also allow the graphene nanoribbon subject to 
two ac fields. In that case, we let j = 1,2 in our general 
formalism in Section |TT] . This case has been studied in 
literature, especially irPSlfor superlattices. Eq.(20) takes 
the form 

oo oo oo 

r—1 ni,n2 — — oo ^i , f2 — — oo 

1 + it(/3 + n x u)x + n 2 L0 2 ) 



from which Eq.(23) follows. We have eliminated the time 



dependence by averaging over the period of the fields to 
find the time-independent current j. In the left hand 
side, we replaced (j(t)) = j, and in the right hand side a 
delta function emerges which ensures that v\ = — ^v 2 . If 
lji = lo 2 , then one must put ax = and a 2 = a so that 
V\ = v 2 . However, we shall generalized this to a case 
of commensurate frequencies. We exemplified the case 
by biharmonic having frequencies which can be periodic 
lu 2 = ulox or non-periodic, io 2 ^ iilo\ with /i = 1,2,.... 
These two cases were studied ir l 17 l 18 l for semiconductor 
superlattices. Defining j 0<r = J|^+T)a ^"=i r£ rsf rs , 



Eq. ( 26 I assumes the form 



oo oo 

1 + it (/3 + [m + fin 2 ]uji) 

Simplifying further, we linearize with respect to one of 
the field amplitudes (say, E 2 ). For a week field, p 2 « 
1 and J n (P) ~ (P/2) 2 /nl, which allows us to take n 2 
(n 2 - v 2 ) = ±1 (0, ±1). We obtained 

OO OO 

r—1 ni— — 00 



Jo(P2) J2n 2 =±l J"2 (p2)Jni (PljJni+lin?. (fi\) 

1 + it(Pq + [m + iin 2 ]ui) 



(28) 



IV. DISCUSSION 

In FigjTJ normalized current density j/jo is plotted 
against reduced static electric field Eo/E cr for aGNR 
(left) and zGNR (right) for an applied ac field. At low 
fields up to Eo(j max ), the quantum derivative of the 
j — E characteristic yields a positive slope. A negative 
slope results for Eq > Eo(j max ). The whole of the region 
E > E (j max ) gives what is called Negative Differential 
Conductivity (NDC). A consequence of NDC in GNRs 
is a formation of electric field domains that impedes a 
continuous motion of electric field waves and thus blocks 
high frequency generation in these nanoribbons. NDC 
disappears quite faster in aGNR as E — ^ 00 as compared 
to zGNR which is more rubust at this limit. 

The curves in Fig(2] demonstrate NDC, they are ob- 
tained at wave mixing of two commensurate frequen- 
cies, uj 2 — fjui. Fig{2] (left) fi = odd and Fig{2] (right) 
/.t = even. The onset of NDC in odd-harmonics oc- 
curs around Eq ~ Ex, and in euen-harmonics it starts 
at Eq < Ex- In both cases, as in the previous NDC 
graphs, lot << 1. 

The combined effect of phase shift and ac amplitude 
on NDC is depicted in Fig(3j There are three peaks at 
low bias fields at points (E ~ E cr , a = 0), (E ~ E cr ,a — 
it) and (E ~ E cri a — 2tt). For now, it is not clear 
what these crests and throughs represents, they might 
be associated with field domains along one direction (for 
a = 0, 2ir) and others along the opposite direction (for 
a — 7r) and vice-versa. 



V. CONCLUSION 

We have demonstrated that graphene nanoribbons ex- 
hibit NDC regions in its j — Eq characteristics at low 
bias field when lot « X. NDC is observed either in 
the presence of bias field alone or by superimposing ac 
field amplitudes on the bias field. For one ac field, NDC 
occurs around lot ~ 0.1. When two ac fields at com- 
mensurate frequencies are applied, high-harmonic NDC 
emerge for both even and odd harmonics at rather very 
low frequencies lot ~ 0.01. The even-series gives pro- 
nounced high-harmonic NDC than the odd-series. The 
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FIG. 1. Negative differential conductivity for (left) armchair and (right) zigzag graphene nanoribbons at different ac field 
amplitudes, ujt = 0.1. The onset of NDC is at low static fields around Eo ~ E. It departs from this condition at high fields 





FIG. 2. Negative differential conductivity due to wave mixing of two ac field amplitudes, (left) fi — 2 and (right) fi — 3. The 
parameters used are E\ = 0.2E cr , E2 = E cr and ujt = 0.01. 



presence of high-harmonic NDC means that it is possi- 
ble for high-frequency generation in graphene nanorib- 
bons when electric field domains are suppressed at high 
enough applied frequencies ujt >> 1 and Eq > E cr . We 
therefore suggest this approach for the study of terahertz 
generation in graphene. 
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